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In designing a computer system, it is vitally important to be able to 
predict the performance of the system. Often, quantities such as 
throughput, processor utilization, and response time can be predicted 
from a closed queueing network model. However, until now the 
computations involved were not feasible for large systems which call 
for models with many processing centers and many jobs distributed 
over many classes. We give a radically new approach for handling 
such large networks — an approach that begins with a representation 
of the quantities of interest as ratios of integrals. These integrals 
contain a large parameter reflecting the size of the network. Next, 
expansions of the integrals in inverse powers of this large parameter 
are derived. For cases in which the number of processing centers is 
greater than one, this is the only technique we know of that yields the 
complete asymptotic expansion. Our method for computing the terms 
of the expansion can be interpreted as decomposing the original 
network into a large number of small "pseudonetworks." Our tech- 
nique also yields easily computed error bounds when only the first 
few terms of the expansion are used. This method has been imple- 
mented in a software package with which we can analyze systems 
larger by several orders of magnitude than was previously possible. 

I. INTRODUCTION 

Closed Markovian queueing networks, which are tractable in having 
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the product form (or separability) in their stationary distribution, 
continue to have a profound influence on computer communication, 
computer systems analysis, and traffic theory. 1 " 4 The closed networks 
have been used to model multiple-resource computer systems, 2,5 
multiprogrammed computer systems, 2,6 " 8 time-sharing, 2 and window 
flow control in computer communication networks; 9,10 networks with 
blocking 1112 require the analysis of a large number of closed networks. 
Not surprisingly, considerable effort has gone into devising efficient 
procedures for computing the partition function, 13 " 17 an element of the 
product form solution requiring significant computation. More re- 
cently, mainly spurred by parallel technological development in com- 
puter communication, there has been a focusing of effort on large 
closed networks 11,1518 " 21 with many classes of jobs and transactions and 
large populations in each class. The point of departure of this effort is 
the realization that the earlier recursive techniques for computing the 
partition function are severely limited in terms of computing time, 
memory storage, and attainable accuracy when it comes to the large 
networks presently demanding analysis. 

In an earlier paper, 2 we introduced a new approach to calculating 
the partition function. We showed there that the partition function 
could be represented as an integral containing a large parameter which, 
in some sense, reflected the large size of the network. In general, the 
partition function is represented by a multiple integral. However, in 
the special case where there is only one node at which queueing can 
occur (a node of type 1, 2, or 4 in the terminology of Ref. 4), the 
partition function is represented by a single integral. In Ref. 21 we 
applied standard theory to obtain asymptotic expansions of the integral 
representation. The standard techniques, however, cannot be extended 
directly to multiple integrals. In this paper, we make use of the special 
properties of our integral representation and obtain a method for 
generating asymptotic expansions of our integrals which works equally 
well for single or multiple integrals. For our single integrals, the 
techniques developed in this paper are easier to apply than the 
standard techniques we used in Ref. 21. For our multiple integrals, our 
technique is the only one we know of which can be used to obtain 
higher order terms in the asymptotic expansion. 

The computational effort of solving a large network with p classes 
and with populations in each class of the order of magnitude of 100 is 
here reduced to be roughly as complex as solving by older techniques 
for the partition functions of (p choose 4) networks where each of 
these networks has a total population of, at most, seven allocated over, 
at most, four classes. Thus, we have reduced the problem to the 
solution of a large number of small problems. One consequence is that 
forp large enough, even our technique will be computationally intrac- 
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table. Our results nevertheless allow large, previously intractable net- 
works to be solved with vastly reduced computational effort and with 
great accuracy. For example, a network with 20 classes, each class 
having a population of 100, can be handled easily by our technique. 
Two noteworthy aspects of this work are (i) the built-in notion of 
depth of accuracy: by computing more terms of the asymptotic expan- 
sion, it is possible to match computational effort to desired degree of 
accuracy, and (ii) a comprehensive error analysis that allows estimates 
to be accompanied by sharp error bounds with little incremental effort. 
Let us elaborate. 

Section III recapitulates and extends the results in Ref. 21 to obtain 
representations as integrals of the partition function of most, but not 
all, closed product-form networks. 4 The class-by-class breakdown of 
the utilization of each processor, itself simply related to mean response 
time and throughput, is given in terms of a ratio of two integrals. 
These are multiple integrals with multiplicity equal to the number of 
centers in the network where queues may form. 

The asymptotic expansions are in powers of (l/N) where N is a 
parameter designed to reflect network size. Our computational expe- 
rience has been that five terms in the mean value expansions are 
generally adequate both for large networks, for which our asymptotic 
technique is particularly well suited, and, to a surprising extent, for 
small networks as well. 21 Section IV gives the procedure for generating 
the general coefficient of the expansion, while the leading coefficients 
are explicitly derived. 

A remarkable feature of the composition of the coefficients make 
their computation amenable to various techniques. As shown in Sec- 
tion 4.4, the coefficients turn out to be very simply related to the 
partition function of a certain hypothetical network which we call the 
pseudonetwork. The topology is related but not identical to that of the 
given network and the processing rates are quite different. Most 
importantly, however, to compute the leading expansion coefficients, 
it is necessary to compute the partition function for only small popu- 
lations in the pseudonetwork. Thus, to compute five terms of the 
utilization expansion, we need only consider the total population over 
all classes to be at most seven in the pseudonetwork. Small population 
in the pseudonetwork has the consequence of its partition function 
being solvable by existing recursive techniques of proven efficacy. 

Section V proves that the series in (l/N) given in Section IV is 
endowed with properties substantially more desirable than those pos- 
sessed by asymptotic expansions in general. It is shown in Section 5.2 
that the truncation error is numerically less than the first neglected 
term of the expansion, and has the same sign. Thus, except for the 
effort in computing an additional term, all that is generally needed for 
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an error analysis is already available in the basic series that is com- 
puted. 

We note here that while he mathematical literature on single 
integrals is extensive, 22 " 24 th re is little on asymptotic expansions of 
multiple integrals. Two notable investigations along lines different 
from here are Bleistein's 25 and Skinner's. 26 Both are marked by extreme 
complexity. Bleistein gives the leading coefficient, while Skinner ob- 
tains the second term. Both terms are quite complicated. 

Let us now elaborate on some limitations of the paper. An important 
conclusion of Ref. 21 is that qualitatively different expansions of the 
integrals exist depending on whether usage is "normal," "high," or 
"very high." This is even more true in the present context of multiple 
processing centers. Therefore, this paper is devoted exclusively to the 
case of normal usage. We propose to consider the remaining usage 
conditions in the future. Exactly what is meant by normal usage is 
explained in Section 4.1. 

The paper also assumes that, for each class of jobs, the routing 
through the network contains at least one infinite server (is) center. It 
turns out that for networks in which this is not true, the asymptotic 
expansions are more appropriately derived in the context of either 
high or very high usage conditions. However, certain basic results on 
the integral representations of partition functions and mean values are 
derived in Section III regardless of whether is centers exist for all 
classes. This paper does not allow load dependent service rates in the 
first-come-first-serve centers. 

The results in this paper have been incorporated in a large software 
package and this will be reported elsewhere. No results on moments of 
individual queue lengths are given here. However, their integral rep- 
resentations and asymptotic expansions are very similar, 21 and the 
details will be published elsewhere. 

II. PRODUCT FORM IN STATIONARY DISTRIBUTIONS: PRELIMINARIES 
2. 1 Product form 

We recapitulate some of the well-known results 4 concerning product 
form in stochastic networks and present them in the form that will be 
used later. 

Let p be the number of classes of jobs and reserve the symbol j for 
indexing class. Hence, when the index for summation or multiplication 
is omitted, it is understood that the missing index is j, where 1 < j < 
p. A total of s service centers are allowed. We will find it natural to 
distinguish the centers of types 1, 2, and 4 which have queueing from 
the remaining centers of type 3 which do not. (The definition of type 
1 to 4 centers is given in Ref. 4.) Thus, centers 1 through q will be the 
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queueing centers, while (q + 1) through s will be the type 3 centers, 
which have also been called think nodes and is nodes. We reserve the 
symbol i for indexing centers. Also, whenever class and center indices 
appear together, the first always refers to class. 

Let the equilibrium probability of finding n,, jobs of class / at center 
h 1 ^ j < p, 1 < i < s, be 7r( yj , y 2 , • • • , y„), where 

yi A {n u , n 2i , • • ■ , n pi ), l<t'<s. (1) 

Closed networks are characterized by conservation of jobs in each 
class. That is, the population of jobs of the y'th class is constant at Kj, 
say. The well-known results on closed networks with the product form 
in its stationary distribution may be given in the following form: 

*Kyi, ••• , y.) =7^n w«(y<)i (2) 



where 



•fo>-G:*)iii(^4i. isisg, 



«/«' 



In the above formulas, we have taken into account the previously 
stated assumption; namely, for the first-come-first-served discipline in 
type 1 centers, the service rate is independent of the number of jobs in 
queue. Also, in (3), 

expected number of visits of class / jobs to center i 

Pji = : 7 7-j . . , . : , (4) 

service rate of class j jobs in center i 

where the numerator is obtained from the given routing matrix by 
solving for the eigenvector corresponding to the eigenvalue at 1. 
In (2), G is the partition function, and it is explicitly 

G(K)= I ... J IjXy,), (5) 

l'n,-/fi Vn —Kj i=i 

where we have written l'n, for £?=i nji and the condition l'n, = Kj to 
indicate the conservation of jobs in each class. Th is, 

g(k)-e-e n (ZMin^ll n (n^jl (6) 

_i=i t «y«!J i-q+i [ nji 1 .] 

In Section III, we will refer to this expression for the partition function. 
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2.2 Asymptotic expansions 
A series 

1 A k /N k 

k=0 

is said to be an asymptotic expansion 22-24 of a function I(N) if 

m-l 

UN) - S A k /N k = 0(AT m ) as N-^ oo (7) 

A-0 

for every m = 1, 2, • • • . We write 

I(N) ~ i A k /N k . 

A=0 

The series may be either convergent or divergent. 

III. INTEGRAL REPRESENTATIONS 

As the representations presented here are basic to the subsequent 
development, we have allowed some duplication in Section 3.1 with 
Section 10 of Reference 21. 

3.1 Partition function 
We start with Euler's integral representation for the factorial, 



n\= e- u u n du. (8) 

Jo 

Returning to (6), we use this representation to write 

(£ nji)\ = [ e~ Ui n u^dm, i - 1, 2, • • - , q. 0) 

Jo 

Substitution in (6) gives 

o-r...f\J-J.) i ... _ [nln^l 



n ng 



dai ••• du q . (10) 



Now by the multinomial theorem, 



g = (n TO" 1 J • ■ • J ex p(- 2 "') 

•II 1 1 ww + £ p;« [ rf-n • • • rfw «7- ( n ) 
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It is noteworthy but not surprising that the parameters p„ for all the 
type 3 centers appear lumped together. Hence, we may simplify the 
notation by introducing pp, where 



P/>= £ Pfi> /'- 1.2, ••• ,p. 



1-9+1 



(12) 



Another consequence of the notation is that the center index i may 
henceforth be understood to range over the processing centers only, 
i.e., 1 < i < q. 

The new quantity pjo has the physical significance of being the 
weighted combination of all the mean think times of the is centers in 
the routing of the y'th class. In particular, if the routing of the yth class 
contains at least one is center, then pjo > and otherwise pjo = 0. Let 
/ be the collection of indices of classes of the former type and let /* be 
the complementary collection, i.e., 



;'E/»p/,>0 and ;'e/*«p >0 = 0. 
With this notation, 



(13) 



Pji 



xnji + E-K.-} n\ipjtui\ dut — dut 

JeJ L i PJO ) jei' I i J 

In vector notation, which we shall use widely, this reduces to 



G = 



UPj6/UKA 
J& / i 



e^IJ^ + ^rfu, 



(14) 



(15) 



- JQ+ 



where* 



u = (Mi, u 2 , • • • , u q )' 

1 = (1, 1, .-- , 1)' 
r J ■ (fji, rfl, • • • , rjg)', 

rji = pji/pjo if JE.I 

= Pj , if ye/* 
Sjj =1 if j e / 
= if jer 

Q + = {u|u>o}. 



i<y<p 



(16) 



• Unfortunately, r„ here is defined to be the reciprocal of the natural extension of r, 
in Ref. 21. 
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We now introduce the large parameter N and define 

PjAKj/N, 1<><P, (17) 

r,A2Vr„ l<y'<p. (18) 

The suggestion in the notation is that in the generic large network 
{fit) and {I;,} are 0(1). That is, the ratio of processing time to think 
time is, in order of magnitude estimation, proportionately less for 
increased populations. There is great latitude in the choice of the large 
parameter. The guiding principle in choosing it should be that the 
resulting values of {/?,} and {Tji) are as uniformly close to 1 as possible. 
In practice, we have used 

JV-max(-l. (19) 



I 1 ! 
tiax < — > . 

ij (rji) 



On substituting (17) and (18) into (15) and after the change of variables 
z = u/N, we obtain from (15) another useful integral representation of 
the partition function which is distinguished by its dependence on N. 
Summarizing for future reference, we have 
Proposition 1: 



G(K) = 



where 







Iwiipft/nK/. 


f e- mu) dz, 
Jo* 


f(z) A l'z - £ 


PA 


Og(8yr + r». 



(20) 
(21) 



(22) 

D 



3.2 Mean values 

We restrict our attention to the mean value K«(K) which gives the 
utilization of the t'th processor by jobs of the ath class for a population 
distribution by class in the network denoted by K = (Ki , K 2 , ••• ,K P ). 
Other interesting mean performance indices, such as throughput and 
mean response time are known to be simply related to {u OI (K)} and 
the interested reader may consult Ref. 17. 

In Ref. 17, 

tf)-p /y , l<o< A l^<7, (23) 

where e„ is our notation for the vector with the oth component unity 
and all other components zero. Thus, the value for the partition 
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function is needed for the given population distribution and also for 
the population in the ath class reduced by 1. 
Now from (20) 



G(K + e„) = 



TT K 

Ml 



e- ru (l + rtu) 



(Kc + DYIK/. JQ . 
X J] (S,y + r/u)*>du if oE.1 



(24) 



0p5 

MI 



-l'u 



«n) 



x [J (5// + r/u)^du if a£/« 



(25) 



From (23) to (25) and the same change of variables, namely z = u/N, 
employed in transforming (20) to (21) we obtain the following repre- 
sentation of the utilization in terms of integrals, 

Proposition 2: For class index a, 1 < a < p, and center index i, 
1 < i < q, 



MK + e.)- 1 - 



1 



r,i(K. + 1) 



Jo* 



Z)e~ mz) dz 



S a i + -2-z . (26) 



,-Nf(*) 



dz 



D 

Some digressionary comments are as follows. Note that the above 
conceals that r OI is normalized differently, but not unexpectedly, de- 
pending on whether o G / or otherwise [see (16)]. Also note that in the 
typical large network for normal operating conditions we always expect 
r„i to be 0(l/iV), precisely because of the normalization used, and K to 
be 0(N) so that the term in braces in (26) is then 0(1). 

IV. ASYMPTOTIC EXPANSIONS 

We henceforth consider only networks in which the route for each 
class always contains an infinite server center. Specifically, 



pjo > 0, j - 1, 2, • 

and the set /* is empty. 

4. 1 The assumption of ' 'normal usage ' ' 
Define 

« 4 1 - 1 m 



,p, 



(27) 



(28) 
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so that in terms of the original network parameters 

«-l-EJ&— , * = 1,2, ...,<?. (29) 

P/o 

It is important to note that a is independent of the choice of N. The 
parameter «,(— oo < a, < 1) is an indicator of the unutilized processing 
capability of the ith center. Positive values of a, correspond to less 
than 100 percent utilization of the processor and negative values which, 
of course, can occur to very high utilizations.* Normal usage in large 
networks will almost certainly require a, > 0, and in all likelihood a, 
will not be close to for all i. 
Hereafter, we assume 

a,>0, i=l, 2, .--,q, (30) 

which condition we refer to as normal usage. Moreover, as a, for some 
i comes close to 0, the expansions given here are not as efficient as 
those derived specifically for such conditions and which we propose to 
give in the future. 

A justification of the usage interpretation that we have given to a is 
provided by a result obtained later (see below Proposition 6) which 
states that, asymptotic with network size, 

m = utilization of ith processor ~ 1 — a, ■ . (31) 

An obvious caveat is that this result is derived for the assumption in 
(30). However, as the utilization in (31) can come close to unity even 
while (30) is satisfied, (31) suggests that for large networks normal 
usage will not extend beyond the range a > 0. 
Observe that 

/•(0) = (32) 

and V/(z) = 1-2—^7-1^ (33) 

so that 

a = V/-(0). (34) 

The assumption of a > and the form in (33) ensures that the function 
/"has no stationary points in Q + since 

V/(z) > V/(0) > 0, z £ Q + . (35) 

Also observe that 



* Unfortunately, a here has an opposite sign from the natural extension of a as 
denned in Ref. 21. 
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^)"^ w (36) 

from which we note that the Hessian is positive semidefinite. 

To conclude, with normal usage, f is a convex function, with its 
minimum in Q + attained at and with no point in Q + where its 
gradient vanishes. 

4.2 Transformations on integrals exploiting normal usage 

Consider the following transformations on the basic integral: 



f e~ Nfw dz - [ 

Jq* a 

-i 



e-" az exp 



-\ NZPjiT'jz -\og(l + rjz)} \dz 



j 



»}]< 



= N~ q e~ au exp 



X/^r;u-Ariog(i + irja 



rfu, (37) 



where u = Nz. Now make the following change of variables, 

Vi A cMii, l<i<q (38) 

and normalize the system parameters with respect to a, thus, 

fji A Tji/m, 1 < j < p, 1 < i < q . (39) 

Observe that in particular 

I> = fjv. (40) 

From (37), 



! 



e- N/lt) dz = -^- e- Vv H(N-\ v)rfv, (41) 

Q+ (Uou) J Q+ 



where 

H(N-\ v) A e °i»-'y\ (42) 

■car 1 , v) a - j /Jfjv - Niog^i + Afyj j. (43) 

Our notation here suggesting N" 1 as the independent variable may be 
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perplexing at this stage, but it provides a clue to the direction of the 
analysis. 

Here are two further digressionary comments: the above transfor- 
mations are meaningful only in the context of normal usage, i.e., a > 
0. In a similar vein, an interpretation of {f 7 ,} as renormalized [Tji) is 
only meaningful for a > 0. It is noteworthy that hereafter we shall be 
dealing exclusively with {Tji} and not at all with {Tji}. 

We need to repeat the transformations given above for the integral 
j(T' z)e~ Nflz) dz. The result may be combined with (42) and (43) to give 
the following compact representation: 

(r a z) m e- mz) dz=^—— e- lv (r o x) m H(N-\v)dv, (44) 

Q+ (n«.-) J Q+ 

m - 0, 1, 2, • • • . □ 

We may now use these expressions in Proposition 2 to obtain for 1 < 
a<p and 1 < i< q, 



u ai (K + e,)- 1 = 



f PoO 1 

\ Pai (K a + l)j 



1 + 



e- lv (T;v)ff(JV-\ v)dv 



N 



I-- 



H(N-\ v)rfy 



(45) 



Note that the bracketed term is independent of the processing center 
index i. 

Let us agree to call 



UN) A I e- lv H(N-\v)dv 



and 



r a l, (N) i 



Jo* 



v)H(N-\ \)dv, 



(46) 



(47) 



where the superscript (1) is a mnemonic for first moment. In this 
notation, we have for future reference 
Proposition 3: 

] (N) 
I(N) 






(48) 

D 



The asymptotic expansions of I(N) and I { „ l) (N) are considered next. 
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4.3 Asymptotic expansions 

This section will outline the procedure for obtaining the asymptotic 
expansions together with plausibility arguments. The proofs of the 
assertions concerning asymptoticity will follow in Section V. Moreover, 
we defer till Section 4.4 certain observations which make feasible the 
efficient computation of the coefficients of the asymptotic expansions. 

Our procedure for I(N) is to first obtain a power series in N~ l of 
H(N~\ v) and then to integrate term by term. Thus, we let 



, ,-, h k (\) 



H(N-\ v)=^ (49) 



and 



and claim that 



A k A e - ly, h k (v)dv (50) 



W)~lj£- (51) 

Let us elaborate on the coefficients {hk(v)) in (49). 

MV) = i W7W m °- V) ' * - °. L 2. • • • ■ (52) 

To make these explicit, we need to first introduce 

A(v) 4 tlL i fi/Jfyrf, k = 1, 2, • . . . (53) 

K j 

Their role becomes clear if we recall (42): 

H(N~\ v) = e s(A, " , - v, > (54) 

and note from (43) that for fixed v E Q + , s(N~\ v) and, hence, 
H(N~\ v) are functions of N~ l , analytic in Re(N~ l ) > e(v), where 
c(v) < 0. Then, 

s <*>(0, v) = -k\f k+1 (y), k - 1, 2, - . . . (55) 

To proceed now to the derivatives of H( • , v) itself, we will find useful 
the following expression in which it is understood that all derivatives 
are with respect to N~ l : 

H^N- 1 , v) - i ( k )s (k+1 - m) (N-\v)H {m \N-\v), 

/ii-o \ m / 

k - 0, 1, • • ■ . (56) 
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From (52), (55) and the above, it is easy to see that there exists the 
following recursive scheme for generating {A*(v)}: 

ho(v) « 1 

1 * 
h h +i(v) = - -r 2 (k + l- m)f k+ 2-m(v)h m (v), 

k + 1 m _o 

• k = 0,1,2,.-. OV/) 

In particular, the leading elements are 
h k (v) - 1, k = 

= -f 2 (v), A = 1 

= -/ 3 (v)+i/1(v), A = 2 

= -/ : 4(v)+/ 2 (v)Mv)-i/i(v), k = 3. (58) 

To summarize the steps discussed so far in the generation of the 
asymptotic expansion, we have 
Proposition 4: 

urn ~ i 4, (59) 

where A* = Je~ rv h k (v)dv, {h k (v)) is obtained recursively from (57) 
with leading elements exhibited in (58), and { /*(v)} is as in (53). □ 

The aspect of the above asymptotic expansion of the integral I(N), 
which consists of decomposing the integrand into the product of an 
exponential and a function, the expansion of the latter function in a 
power series, and the final term-by-term integration is like the proce- 
dure which, in the context of single integrals, is justified by Watson's 
Lemma 23 under certain conditions. Our contribution has been to show 
that a generalization of this fundamental result exists for the multiple 
integrals of interest in stochastic networks. 

Our procedure for obtaining the asymptotic expansion for Io ] (N) is 
very similar and consists of obtaining a power series in N" 1 of 
(t' y)H(N~\y) and integrating term by term. However, we notice the 
simplifying fact that 

(f..v)ff(N-,v).E^^. (60) 

k-0 M 

Thus, the procedure for this integral is as follows: 
Proposition 5: 
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where 






00 A {1) 

A=0 ^V 



e- iv (r;v)^*(v)dv, * - 0, 1, 2, ... . 



(61) 



v* 



(62) 

D 
As observed in Ref. 21, the asymptotic expansion for ..he integrals may 
be used to generate asymptotic expansions for their ratios on account 
of powers of N' 1 forming a multiplicative sequence. 22 Thus, the coef- 
ficients of an asymptotic expansion for u OI (K + e„) _1 may be obtained 
from formal substitution in Proposition 2 of the expansions in Propo- 
sitions 4 and 5. This gives 

Proposition 6: 



Poi(K a + 1) 



PoO 



(1) 



u ai (K + e o y 






A 



1 

A~ 



k = 



in 



A a ,k 2j A m B a ,k-m 



k - 1, 2, 



(63) 



(64) 



□ 



With the above proposition we may generate (£ + 1) terms of the 
expansion for «„, from k terms of I(N) and /^(iV). 
An immediate corollary to the above proposition is 

PoiK a 



u oi (K) ~ 



PaO 



and 



«,(K) = utilization of ith processor = £ w y ,(K) 



~l-a„ i = l,2, .-., 9 , (65) 

which was claimed earlier in Section 4.1 in the course of giving physical 
meaning to the parameters {a,}. 

This corollary illustrates the important point that the terms in the 
asymptotic expansions A k /N k , A l a ]l/N k , and B aJl /N k+l are all inde- 
pendent of N and depend only on the network parameters. The dummy 
parameter N serves to show how to group the terms of the same 
magnitude. This independence from N follows from the fact that a is 
independent of N, as noted earlier, and from the fact fk(v)/N k ~ l is 
easily seen to be independent of N. Once this is noted, the result 
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follows easily from the definitions of the terms. The choice of N is 
important numerically, which will be discussed in a subsequent paper. 

4.4 Pseudonetworks and the computation of expansion coefficients 

Here, we consider the compositions of coefficients, {A k } and {A ok }, 
and find that quite remarkably they are related intimately to the 
partition function of a certain hypothetical network which we call the 
pseudonetwork. It turns out that to compute the leading elements of 
{Am} and {A { Jj,}, we need to consider the pseudonetwork with only 
small populations. Thus, existing techniques known to be effective for 
computing partition functions for small populations may be used to 
compute the leading coefficients of our asymptotic expansions. 

An example will prove useful. From Proposition 4 and eq. (58) we 
see that 

A 3 = J e- x A -/ 4 (v) + /- 2 (v)/3(v) - i/1(v) Idv. (66) 

Now consider only the third term after denoting it by A 3 a. Using the 
expression for f 2 (v) as given by (53), we obtain 

A 33 = -^Z& [ e" iv (f;v) 6 rfv 

-■hi hlPh [ e- lv (f>) 4 (f ' k v) 2 dv 

lo Mk J 

-ill A&& I e- 1,v (r;v) 2 (nv) 2 (^v) 2 rfv, (67) 
48 j+k+i+j J 

where the subscripts j, k, and / are class indices ranging over [1, p]. We 
now make the observation that the generic integral in the composition 
of the asymptotic expansion coefficients is within a multiplicative 
constant of 



g(m) = g(m h 77i2, • • • , m p ) A 

'n 



— v f e- lu n(f>Pc?u. (68) 

9-m) Jq+ J 



The above is an important form for we may now identify it with 
quantities previously encountered. 

We first give the following equivalent expression for#(m). 

l'ni-mi l'np-mp i-1 I \ j / j 'f/«- J 

Recall the expression in Proposition 1, eq. (20), for the integral repre- 
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sentation of the partition function. Specialize the expression there to 
a network with no infinite server centers, i.e., set / is empty, and find 
that it reduces to 



G(K)=- e- lu IT(r;u)*tfu. (69) 



On comparing (68) and (69), or (68') and (6), we may conclude that 
g(m) is the partition function of a certain network. 

Call this hypothetical network the pseudonetwork. What character- 
izes the pseudonetwork? To begin with, it is closed and lacks is centers. 
There are, as in the original network, exactly q processing centers and 
p classes of jobs. The processing rate of jobs from the y'th class in the 
jth center of the pseudonetwork is Iy, , where, you will recall, fy, = 
Tji/cti. In agreement with past convention, (mi, m 2 , • • • , m p )' denotes 
in vector form the population distribution by class in the network. 

We may follow the procedure outlined in the example concerning 
A33 above to express all the leading coefficients Ak, k ■• 0, 1, 2, 3 in 
terms of the partition function of the pseudonetwork. This gives 

A = 1 

j 
A 2 = 2l frgttej) + 3 J $fg(4ej) + i J fofagitoj + 2e k ) 

J j ^ Mk 

A3 = -6 2 0jg(4ej) - 20 2 fi*g(5ej) - 15 £ tfgitej) 

i J j 

-2 2 fijPkg(2ej + Se h ) 
- 3 2 tfPkgdej + 2e k ) 

-« 2 foPkPig(2ej + 2e k + 2ei). m 

In these expressions, j, k, and I are class indices each with range 
[1, Pi 

Notice that in the computation of A k , k = 0, 1, 2, 3, the population 
distribution m that appears in #(m) may be characterized as being 
quite small. The total population in the pseudonetwork over all classes, 
2 rrij, is at most 6. A further simplifying condition is that we need 
consider only distributions with three classes at most with nonzero 
populations — the extreme distribution arises in #(2e y - + 2e* + 2e/). 

Since we are interested in population distributions where all but a 
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small number of classes have no members at all, we may equivalently 
choose to view the pseudonetwork as a collection of smaller networks 
each with a small (less than p) number of classes of jobs. For example, 
g(2ej + 3e A ) may be viewed either as the partition function of the 
pseudonetwork in which all but classes j and k have zero population, 
or as originating from a network with only two classes with population 
distribution (2, 3) but one which is otherwise unchanged. Such distinc- 
tions, while not material to the procedures given here, may be conse- 
quential in the efficiency of the computations. 

The coefficients {A { ]1} of the expansion for tfP(N) may also be 
expressed in terms of the partition function of the pseudonetwork. As 
the derivation is similar, it will suffice to give the results, which we do 

in the Appendix. 

It will be observed that to compute {A£*} , k - 0, 1, 2, 3 we need to 
consider various allocations to classes of a total population in the 
pseudonetwork of at most 7. Thus, to compute the leading five terms 
of the utilization u oi , we need A k and {A$} for k = 0, 1, 2, 3, and these 
are computed from the values of the partition function of the pseudo- 
network for various allocations to classes of a total network population 
of at most 7. This is an elaboration of a claim made in the Introduction. 

V. ERROR ANALYSIS AND PROOF OF ASYMPTOTICITY 

Equations (59) and (61) contain claims requiring proof on the 
asymptoticity of the expansions of I(N) and J), 1 * (N). This is provided 
here as a corollary to a complete error analysis. In fact we show that 
the expansions given earlier have properties more attractive than that 
required of asymptotic expansions. For instance, asymptoticity re- 
quires that errors incurred in the estimation of the integrals from the 
use of, say, m leading terms is of the same order as the (m + l)th term 
as n —> oo. We prove that the expansions derived have the stronger 
property that the truncation error is bounded by the (m + l)th term. 
The practical benefit of this analysis is that with very little incremental 
effort we can accompany our estimates of the mean values with sharp 
estimates of the estimation error. 

5.1 Completely monotonic functions 2A 

We need the following definition: for any nonnegative R, let Q + (R) 
be the set of nonnegative vectors with norm bounded by R, i.e., 
Q + (R) = {v|v ^0 and ||v||<i?}. Note that Q + (i?)->Q + as #->oo. 

The proposition below (cf. Ref. 21) states a remarkable property of 
the function H(N~\ v) which is a key to much of the error analysis. 

Proposition 7: For all v G Q + (R), R<*>, H(N~\ v) is a completely 
monotonic (or alternating) function of AT" 1 . That is, 
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(_1) dO/W { ' and 0<ATi<oo. < 71) 

Proo/: Consider the form for i/(AT\ v) from (42 to 43), Section 4.2: 

# (AT 1 , v) = J] | e"^ v) (l + i Pjv\ . (72) 

Because products of completely monotonic functions are also com- 
pletely monotonic, it suffices to show that 

(l + ilfrr)*- 

is a completely monotonic function of N~\ Let us write 

(l+if/vW = e' (w> , (73) 

where 

t(w)- — \og(l + aw), (74) 

w 

by identifying w = 1/N and a = Iyv. We note that < w < oo and that 
a is nonnegative and bounded, with the latter property being ensured 
by the restriction of v to Q + (R), R < oo. Thus, all derivatives of 
log(l + aw), and consequently of t(w), exist and are continuous for 
< w < oo. Since 






rfa;' 



{e MH "} 



(75) 



we may conclude from a simple inductive argument that 

if t(w) is a completely monotonic function of w, 

then e"'" 1 is a completely monotonic function of w. (76) 

Finally, to show that t(w) is a completely monotonic function of w is 
to show that (log(l + x)}/x is completely monotonic. This is true, but 
we omit the proof. □ 

We need the following analogous property in connection with the 
integrand of /^(N): 

Forall v G Q + (R), R< oo, (f 'v)/f (N _1 ,v) is completely 1 

monotonic in N~\ J (77) 

The proof is immediate from the preceding proposition since the 
additional factor (F„v) does not depend on N~\ 
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5.2 Error bounds 

Proposition 8: For all N, < N < oo 

A m ~ i Ah 

^ <m) ~SoW <0 - m = 1 ' 3 ' 5 ---- 

m ~ 1 A L A 
0</(N)-S^<^ m = 2,4,6, .... (78) 

Proof: We initially require v E. Q + (R), R < ™ so that the preceding 
proposition is applicable. Viewing H{N~\v) as a function of AT 1 , we 
may use a version of Taylor's theorem 24 that is accompanied by an 
estimate of the truncation error for the series to obtain 

w^'Z^+wizmNr™^ (79) 

where 

|e[0,N- ] ]. 

Consider first the case of m odd. From the preceding proposition, the 
mth derivative of H(N~\ v) is a nonpositive and monotonically 
nondecreasing function of N~ l . Hence, 

^(0,v)<— — — H(|,v)<0, £>0. (80) 



d(l/N) m ' d(l/N) 

Substituting in (79), 



^ mN -, v) -%^0. (81, 

Hence, 

< ' 7/(V '.vldv- V^* I e- lv A*(v)rfv 

<0 

h m (v)dv. (82) 



(?+(*) * =0 " -V<*> 






N 

The pair of bounds holds uniformly in R. Consequently, we may let 
R — » oo and drop the restriction on R to obtain (78). 

The proof for m even is very similar with the starting point being 
the following replacement for (80), 

•j/n a"i 

^WlNr Hi ^ )s WlNr mo - vh (83) 

□ 
The rest of the proof is omitted. 
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The above proposition states that the error incurred from using only 
a certain number of leading terms of the expansion for I(N) is numer- 
ically less than the first neglected term of the series and has the same 
sign. 

Another implication that can be quite useful in practice is that the 
estimate with an odd number of terms is an upper bound on the 
integral and an even number of terms gives a lower bound. Thus, the 
error sequence alternates in sign. (It is also true but less consequential 
that the terms of the expansion also alternate in sign.) 

In particular the above proposition proves the asymptoticity of the 
expansion in Section 4.3. 

By a matching argument and with recourse to the complete mono- 
tonicity of (f!,v)H(N~\ v), see (77), we also have 

Proposition 9: For any class index a and all N, < N < oo, 
±£ S I»\N)-Z±*X0, m- 1,3,5,... 

OsI^m-l-^S^, m = 2,4,6,.... (84 ) 

With error estimates available for both I(N) and I { a 1] {N), it is 
straightforward to use these to obtain an error estimate for the mean 
value given in Proposition 3, eq. (48). 

APPENDIX A 

The Coefficients {A $ } in Terms of the Partition Function of the Pseudo- 
Network 

Here a is a given fixed class index, while j, k, and / are also class 
indices each ranging over [1, p]. It is also understood that j, k, I, and 
a are all distinct. 

A { l l = g(e ) 

A$ = -3&s(3e ) - I jg (e o + 2e,) 

A3 = 8/? £(4e ) + I5(3ig(5e a ) 

+ I [2/3 jg (3ej + e„) + 3##(4e, + e ) + 3p o p jg (3e + 2e,)] 

j 

+ -'ZPjPkg(2e J + 2e k + e a ). 

A$= -30p a g(5e n ) - 120/S 2 o g(6e n ) - 10501 g (7 e a ) 

- £ [6&£(4e, + e ) + 20/3jg(e o + 5e y ) + 6/3 a /3jg(3e + 3e,) 
/ 

+ 2p o /3jg(4e a + 2e,) 
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+ I5fi)g(e + toy) + 15filfijg(Ge. + 2e,) + dfijjgito. + 4e,)] 



2PjPkg(e + 2e, + 3e*) + 3^j/3 k g(e + 4e y + 2e k ) 



-I 
M L 



+ 5 faPjPkg(3e + 2e> + 2e*) 

- s S PjPkPig(2ej + 2e k + 2e, + e„). 

6 /.A,/ 
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